Monodisperse self-assembly in a model with protein-like interactions 
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We study the self-assembly behaviour of patchy particles with 'protein-like' interactions that can 
be considered as a minimal model for the assembly of viral capsids and other shell-like protein com- 
plexes. We thoroughly explore the thermodynamics and dynamics of self assembly as a function of 
the parameters of the model and find robust assembly of all target structures considered. Optimal 
assembly occurs in the region of parameter space where a free energy barrier regulates the rate of 
nucleation, thus preventing the premature exhaustion of the supply of monomers that can lead to 
the formation of incomplete shells. The interactions also need to be specific enough to prevent the 
assembly of malformed shells, but whilst maintaining kinetic accessibility. Free-energy landscapes 
computed for our model have a funnel-like topography guiding the system to form the target struc- 
ture, and show that the torsional component of the interparticle interactions prevents the formation 
of disordered aggregates that would otherwise act as kinetic traps. 

PACS numbers: 81.16.Dn,87.15.ak,87.15.km,81.16.Fg 



I. INTRODUCTION 

The self-assembly of simple building blocks into larger, 
ordered structures is a ubiquitous process in biology, and 
holds great promise for applications in materials science 
and nanotechnologyP^ In particular, proteins are the 
building blocks for a vast array of biological structures, 
including capsules, fibres, tubes, sheets, channels and 
motors. Here, we restrict our interest to those exam- 
ples where the self-assembly is both one-component and 
monodisperse, i.e. the structures formed are of a specific 
size and made up of multiple copies of the same pro- 
tein. The archetypal examples of this type are icosahe- 
dral virus capsids, which are designed to encapsulate the 
genomic material of the virus. The capsid proteins are 
produced in large quantities inside host cells, and need 
to assemble correctly and spontaneously on a biological 
time scale in order for the virus to propagate. While 
many capsid assembly processes rely on interactions be- 
tween the capsid p roteins and scaffolding proteins or nu- 
cleic acids I214151S1ZI capsids for a number of viruses have 
been sho wn to ass emble in vitro in the absence of these 
moleculesPininini 

There are also a significant number of non-viral pro- 
teins that form sh ell-like homomeric complexes with high 
symmetry! 13 l 14 l 15 l For example, ferritin is a complex that 
is made up of 24 sub-units with octahedral symmetry 
that stores iron inside its central cavity,^ and dihy- 
drolipoyl acetyltransferase can form a 60-unit dodeca- 
hedral complex that is at the core of the pyruvate de- 
hydrogenase multienzyme complex. 17 However, the self- 
assembly behaviour of such complexes has been much less 
studied than for virus capsids. 

Clearly, establishing a good understanding of the na- 
ture of the self-assembly process of protein complexes 
and virus capsids, as well as being of fundamental inter- 



est, will be of value to many potential biomedical and 
bionanotechnological applications, e.g. in the des ign of 
drugs which interfere with capsid assembly 18 * 19 * 20 * or the 
use of capsids as vehicles for the delivery of drugs or 
other agents PUSH Finally, an understanding of the fea- 
tures that enable such robust assembly in these protein 
systems may be very valuable in designing synthetic self- 
assembling systems r 3 ' 24 ' 

There is now a considerable body of work 
studying the assembly of virus capsids, both 

experimenta ll y ' 7 ' 11 ' 12 ' 25 ' 26 ' 27 ' 28 ' 29 ' 30 ' 31 ' and 

theoretically ! 3 ^ 

The transient nature of the intermediates present a par- 
ticular challenge to the experiments. By contrast, 
characterizing species with short lift-times is much less 
of an obstacle for simulation and modelling, and so 
they can play an important and complementary role 
by illuminating the mechanisms of assembly. Instead, 
the problem for simulations is the wide range of time 
and length scales associated with capsid assembly. For 
example, although all- atom simulations of small viruses 
are now possible, 53 the time scales associated with 
self-assembly are far beyond what is currently feasible. 
Therefore, it is necessary to use much simpler coarse- 
grained models, the hope being that if self assembly has 
general rules as to what conditions and subunit designs 
lead to successful assembly, these should be accessible 
through such models. 

In this vein, a number of computational approaches 
have been directed towards capsid assembly. Kinetic 
models consider the populations of particular clusters 
and model cluster growth by assigning rate constants to 
combination/breakup events, using eithe r a d ifferential 
equation^ 3 - 6 - 1 or a discrete event approach JmEU However, 
these models do not include any information about the 
spatial positions of subunits, and only consider a lim- 
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ited set of assembly pathways. Simulations do not suf- 
fer from these restrictions, but it is more challenging to 
generate good statistics and thoroughly analyse the pa- 
rameter space. Indeed, it is only relatively recently that 
models with fully reversible dynamics have been stud- 
ied. The simulations can be grouped into two types de- 
pending on whether the protein subunits are modelled 
as anisotropic objects con sisting of a rigid structure of 
smaller spheres^- 4 * * 46 * 47 * ^ or as spherical particles with 
anisotropic 'patchy' interactions.^ The latter approach 
is typically computationally less intense, and is the one 
we use here. 



II. METHODS 



Model 



The model system we u se is a s lightly modified ver- 
sion of that used previousl y 54 * 60 * 61 * and in the accompa- 
nying paper. 55 It has also recently been used to study the 
assembly of tetrameric protein complexes.^ The model 
consists of spherical particles with a number of sticky 
patches, which are defined by patch vectors. The inter- 
action potential is pairwise and is based on the Lennard- 
Jones potential, 



Previously, 54 we have used such a patchy particle 
model to study the assembly of icosahedral clusters, and 
in the accompanying paper this study is extended to 
the remainder of the Platonic solids. 55 However, in these 
studies the model included no constraints on the torsional 
(dihedral) angle between interacting particles. While 
this is likely to be appropriate for synthetic 'patchy' col- 
loids and nanopart icles JMMl ^ does not represent protein- 
protein interactions well since the complex structure of 
the interfaces between the proteins will tend to restrict 
rotation. Such torsional constraints are also relevant to 
recently created 3- and 5-armed DNA building blocks 
that can assemble into tetrahedra, dodecahedra, trun- 
cated icosahedraj^ cubes^ and icosahedraj^ and where 
the relative orientation of these units is controlled by the 
number of turns of the DNA double helix along each arm. 



Here we consider the effects of adding a torsional 
component to the interparticle potential of this minimal 
model, and proceed to map out the behaviour of our sys- 
tems over a wide region of parameter space. We examine 
the dynamics and mechanisms of assembly, and consider 
the behaviour of the system in regions in which assembly 
is not successful in order to understand those processes 
which compete with successful assembly. This work has 
a number of distinctives compared to previous simulation 
studies of viral assembly. Firstly, we do not just consider 
the assembly of an icosahedral target that is equivalent 
to a T = 1 capsid, but more generally consider structures 
with other symmetries that are also relevant to non- viral 
protein complexes. Secondly, the model is simpler than 
those previously studied, thus allowing us to survey the 
behaviour of the model very comprehensively. In particu- 
lar, we are able to connect the dynamics to the underlying 
thermodynamics of the model, and obtain detailed pic- 
tures of the free energy landscapes governing assembly. 
Finally, comparison to the non-torsional model consid- 
ered previously, allows us to understand what features 
of the behaviour arise because of this component in the 
interparticle interactions, and such information will be 
particularly relevant for those trying to design synthetic 
particles that are able to self- assemble. 
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but the attraction is modulated by a pair of orientation- 
ally dependent terms, V ang and V tor . V ang is a function 
of how directly the two patches are pointing at one an- 
other, while the additional factor V tor is a function of the 
torsional angle between the two particles, i.e. it varies as 
one of the particles is rotated about the vector connecting 
the two particles. Thus, the complete potential is 

v l3 (r l3 ,n u n 3 ) = (2) 

f V L j(rij) r < cr LJ 

\ VLjiri^V^irij.ni.n^Vtoririj^i^j) r>a L j, 

where is the orientation of particle i. For any pair of 
particles only the pair of patches that maximizes V ang V tor 
are considered to interact. V ang has the form: 

^ang(?ij, fti, flj) = Gijftij, Qi)Gji(rji, Qj) (3) 



Gij(Yij,Q,i) 



exp 



^ang 
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where <r ang is a measure of the width of the Gaussian, 
Okij is the angle between patch vector k on particle i and 
the interparticle vector r^-. 

In the calculation of Vt or , an additional reference vector 
on each particle needs to be defined. We have chosen this 
reference vector in each case to lie on the symmetry axis 
of the particle. V tOY is a maximum when the projections 
of the two reference vectors onto the plane perpendicular 
to the interparticle vector lie parallel. Specifically, if <j) is 
the angle between the projections, then 



exp 



2<x t 2 or 



(5) 



where a tOY is a measure of the width of the Gaussian. 
The effect of the inclusion of Vtor in the potential is to 
penalise twisting around the interparticle vector, with 
smaller values of <T t0 r giving a stronger constraint. In 
order to reduce the number of parameters to consider and 
because we consider it physically reasonable that cr t0 r and 
<j ang are coupled, we set <7 t0 r = 2<j ang throughout this 
paper, so that the specificity of interactions is given only 
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in terms of <7 ang . Our results are relatively insensitive to 
the precise value of the ratio of these parameters. 

The patch vectors are chosen such that they point di- 
rectly at the neighbouring particles in the target struc- 
ture. This choice, along with that of the reference vector 
above, ensures that the target structures are lowest in 
energy for sufficiently specific patches. The target clus- 
ters that we consider are the five Platonic solids. Thus, 
each particle has internal C n symmetry because of the 
presence of n equivalent and regularly- arranged patches. 
Therefore, when comparing to virus capsids and other 
large homomeric protein complexes, the particles do not 
represent individual proteins, which have no symmetry, 
but cyclic protein complexes that can further assemble 
into larger complexes. For example, the icosahedron- 
forming particles could be considered to represent the 
pentameric capsomers often posited as intermediates in 
viral capsid assembl y 26 * 29 * and the complete icosahedron 
a T = 1 capsid. Similarly, the dodecahedron-forming 
particles could represent dihydrolipoyl acetyltransferase 
trimers.^ In this context, we are modelling the second 
stage of a hierarchical self-assembly process.^ 

For computational efficiency we use a cut-and-shifted 
version of the potential with a cutoff distance of 3(Jlj, 
and also shift the crossover distance in Eq. [2] so that it 
occurs when the cut-and-shifted potential passes through 



B. Simulation 

All our simulations are based on Metropolis Monte 
Carlo simulations in the canonical ensemble. We restrict 
the translational and rotational moves to be local, so that 
the particle motion mimics the diffusive behaviour of pro- 
teins and nanoparticles in solution. Where equilibrium 
statistics are needed we make use of the umbrella sam- 
pling technique, applying an order parameter dependent 
bias to facilitate the crossing of free energy barriers. For 
more details see the accompanying paperP^ 

A particle number density of 0.15 a^j is used for all 
simulations. This represents a significantly higher con- 
centration than would normally be considered in in vitro 
studies involving proteins. Because the time scales acces- 
sible to simulation are far shorter than those available to 
experiment, it is necessary to use higher concentrations 
in order to observe assembly. For our original model 
without torsional constraints we found that for concen- 
trations in the range 0.025 — 0.4 a^j the final yields only 
show a weak dependence of on concentration.^ 



III. RESULTS 

In order to map the behaviour of our systems over a 
wide region of parameter space, we performed large ar- 
rays of simulations with varying values of the patch width 
<j ang and the temperature T, for each of the Platonic 
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FIG. 1: (Colour Online) The percentage yield of target clus- 
ters formed after 80 000 MC cycles as a function of the patch 
width dang (measured in radians) and temperature at a num- 
ber density of 0.15 cr^j for systems of 1200 particles designed 
to form (a) tetrahedra, (b) octahedra, (c) cubes, (d) icosahe- 
dra and (e) dodecahedra. The insets show the equivalent plots 
for simulations without torsional constraints, where the axes 
span the same parameter ranges as the main plots. No inset is 
included for dodecahedra since without torsional constraints 
there was no conditions under which dodecahedra assembled. 
The images in the top right of each plot show the relevant 
target structure. 
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solids. The results are shown in Figs. [T] and |2j Note that, 
although we will typically talk about the dependence of 
the behaviour on temperature, we could equivalent have 
used the interaction strength, where decreasing temper- 
ature corresponds to increasing interaction strength. In 
Fig. [T] which gives the yields of successfully assembled 
clusters, a region of high yield is visible for each target 
structure at moderate values of temperature and patch 
width, where the target is thermodynamically stable and 
kinetically accessible. Moving away from this region, 
yields fall away as competing processes become domi- 
nant. 

The high temperature region corresponds to a gas of 
monomers and small clusters. Above a certain temper- 
ature, which we term T c , the target clusters cease to be 
stable with respect to this high entropy gas. This tem- 
perature is a strong function of the patch width cr ang since 
wider patches lead to higher entropies for the target clus- 
ters, arising from internal vibrations. Fig. [2] which shows 
the average cluster sizes achieved as a function of patch 
width and temperature, clearly shows the extent of the 
monomer gas region. 

In the early stages of the simulations below T c , the 
particles tend to rapidly assemble into curved shell-like 
structures. The inclusion of torsional constraints in the 
potential favours a certain curvature in the structures 
formed, such that even erroneous structures are typically 
still hollow shells of roughly the correct size and shape. 
In particular the constraints tend to enforce convexity on 
the structures, greatly restricting the range of possible 
structures which can compete with the target structure. 

At high values of the patch width <j ang , the interactions 
between particles become less specific. Low energies can 
still be obtained by significantly distorted clusters, and 
so with larger <j ang increasingly malformed clusters are 
observed, until correctly formed clusters become a rar- 
ity. These malformed clusters are somewhat similar to 
the monster particles which have been observed experi- 
mentally in Turnip Crinkle Virus pa and in simulations 
by Nguyen et al However, the nature of our interac- 
tion potential limits the range of structures which can be 
formed, so that the fused-shell structures seen in these 
papers are not seen. Instead we find convex structures 
similar to those found in more recent work by Nguyen 
and Brooks.^ At still higher values of a ang the patches 
are sufficiently wide that they lose their specificity, and 
instead of forming hollow shells, the system forms large, 
roughly spherical liquid-like droplets. 

At low temperatures more complicated behaviour is 
observed, with generally lower yields but some success- 
ful assembly even at very low temperatures for moderate 
values of a ang . The low yields arise because rapid nucle- 
ation leads to the removal of almost all the monomers 
from the system by the growing clusters, leaving many 
clusters half- formed and unable to grow further. This 
phenomenon of monomer starvation is a recurring theme 
of virus assembly, and has been reported in a number 
of previous experimental studies as well as being pre- 
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FIG. 2: (Colour Online) The mean cluster size (averaged over 
particles) of systems designed to form (a) tetrahedra, (b) oc- 
tahedra, (c) cubes, (d) icosahedra and (e) dodecahedra, for 
the same simulations as Fig. [I] The insets show the equiva- 
lent plots for simulations without torsional constraints, where 
the axes span the same parameter ranges as the main plots. 






FIG. 3: (Colour Online) Snapshots of typical configurations of 72 icosahedron- forming particles, (a) T = 0.08, a ang = 0.35, 
showing monomer starvation, (b) T = 0.08, cr an g = 0.6, showing misformed clusters, and (c) T = 0.16,cr an g = 0.75, showing 
liquid droplets. 



dieted by kinetic modelePEH and simulations! 34 ' 35 ' 43 1 46 1 
We shall return later to consider the reason for the re- 
covery in yield at moderate patch widths. 

Fig. [3] shows typical configurations for small systems in 
different regions of parameter space for the icosahedron- 
forming particles. Fig. |3ja) shows a system displaying 
monomer starvation. Six partially formed clusters are 
visible, as well as one trimer and one 14-particle clus- 
ter which is not easily able to rearrange because of the 
low temperature and patch width. This system is un- 
likely to form any correctly formed clusters within ac- 
cessible timescales. Fig. j3^b) shows a system at higher 
patch widths, such that misformed clusters are relatively 
low in energy and are frequently observed. Two cor- 
rectly formed icosahedra are visible, along with a num- 
ber of misformed clusters. Finally, Fig. [3^c) shows liquid 
droplets formed at very high cr ang . Unlike the misformed 
clusters seen in Fig. [3^b), these droplets are not hollow, 
and the particles are mobile within the droplets. 

The general form of the behaviour of the systems is 
to first order independent of the target structure. The 
plots for each of the targets in Figs, [I] and [2] have es- 
sentially the same form in each case. However, some 
differences are also observed. Most strikingly, the tem- 
perature delineating the transition from target clusters 
to a monomer gas, T c , varies considerably. In general T c 
is greater for smaller clusters (since there is a lower en- 
tropic cost for the formation of these clusters), and also 
greater for higher numbers of patches per particle (since 
the clusters are more energetically favoured). These ef- 
fects lead to fairly similar values of T c for the tetrahedron, 
octahedron and icosahedron, where the differences from 
the two effects largely cancel, and lower values for the 
cube and the dodecahedron. 

The insets in Figs, [I] and [2] show the results for equiv- 
alent systems without torsional constraints. The most 
striking differences are visible in Fig. [2j where the in- 
sets show that in the absence of torsional constraints, 
extremely large clusters are formed over wide ranges of 
parameter space. These large clusters correspond to dis- 



ordered kinetic (at low T) or thermodynamic (at high 
c"an g ) aggregates. Torsional constraints prevent the for- 
mation of these large and disordered aggregates by en- 
forcing convexity, the only exception being at the largest 
values of cr ang in Fig. [2] that are well away from the region 
of successful assembly. 

For systems without torsional constraints, the compe- 
tition between aggregation and self-assembly can severely 
limit the yield of the target structured The effects of this 
competition are particularly noticeable for cubes, where 
assembly is limited to a small region, and dodecahedra, 
where assembly does not occur at all. For these two tar- 
gets the yields dramatically increase when torsional con- 
straints are included because the competing aggregates 
cannot form. 

For the targets that readily assemble with and with- 
out torsional constraints, there is one region in which 
the yield is clearly decreased by the inclusion of tor- 
sional constraints. At patch widths slightly higher and 
temperatures slightly lower than the optimum, the sys- 
tems without torsional constraints assemble by a "bud- 
ding mechanism" P=l In this mechanism the particles first 
coalesce into disordered aggregates, which then rearrange 
and bud off completed clusters. Since torsional con- 
straints prevent the formation of large aggregates, they 
have the effect of inhibiting this indirect assembly mech- 
anism. 

These effects are clearly visible in the different shapes 
of the regions of successful assembly with and without 
torsional constraints in Fig. [I] The "lobe" of success- 
ful assembly due to budding is missing in the simulations 
with torsional constraints, while a new region of partially 
successful assembly arises at low temperature and mod- 
erate patch widths, in part due to the lack of competition 
with large kinetic aggregates. 

Assembly in systems with torsional constraints pro- 
ceeds exclusively by the stepwise addition of monomers 
and small clusters onto growing clusters. Fig.[4]shows the 
average cluster size as a function of time at different tem- 
peratures, with and without torsional constraints. While 
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FIG. 4: (Colour Online) The mean cluster size (sampled over 
particles) as a function of number of MC cycles at a ang = 0.45 
and at different temperatures (as labeled) for a system of 1200 
icosahedron-forming particles. The dashed lines show the 
equivalent results for systems without torsional constraints. 
Each line is an average over ten simulations. 



in the absence of torsional constraints the average cluster 
size sometimes approaches the target size from above (in- 
dicating that assembly proceeds by the budding mecha- 
nism) , the average sizes for the simulations which include 
torsional constraints all approach the target size from be- 
low, consistent with a stepwise growth mechanism. 

Figure [5] provides a detailed picture of the dynam- 
ics by showing the numbers of particles in clusters of 
a given size as a function of time, under four different 
sets of conditions for the example of the icosahedron. 
Under optimal assembly conditions (Fig. |5ja)), there is 
only a small population of intermediate-sized clusters and 
as time progresses the population of monomers falls off 
slowly as more complete icosahedra are formed. That the 
population of intermediate-sized clusters is always less 
than monomers or icosahedra (after an initial lag time 
required for icosahedra to start forming) is indicative of 
a significant free energy barrier for the formatio n of the 
target structure. As has been noted previouslyj 6 | 7 | 32 | 41 l 
relatively slow nucleation is a prerequisite for successful 
capsid assembly. When few nuclei are formed, the popu- 
lation of monomers remains high for longer allowing the 
majority of the nuclei to grow to completion. 

At lower temperatures, such as in Fig.^b), the effect 
of monomer starvation is evident. Under these condi- 
tions, the formation and growth of cluster is all down- 
hill in free energy. The consequent rapid nucleation of 
clusters leads to a rapid decrease in the population of 
monomers, and many of the partially formed clusters are 
unable to grow to completion. At this point assembly is 
effectively stalled, and can only proceed through a slow 
process of breakup of the existing clusters. This is the 
reason for the low yields observed at low values of T and 

^"ang- 

For these combination of reasons, optimal assembly oc- 



curs fairly close to T c . This result is consistent with Zlot- 
nick's assertion that "weak" bonds are sufficient and in- 
deed optimal for efficient assembly of viruse^ 3 - 7 -^ (but of 
course while still being strong enough that the capsids 
are stable) 

Although the phenomenon of monomer starvation is 
found at low temperature regardless of a ang , yields are 
found to recover for moderate patch widths. This recov- 
ery is observed in Fig. [5|c) and is due to an interesting 
mechanism of combination and rearrangement of clus- 
ters. Pairs of partly formed clusters occasionally collide 
and stick together. Most often the two parts will not 
fit perfectly to form a complete cluster, and so a pro- 
cess of rearrangement then occurs. If excess particles are 
present, they may often be ejected from the rearranging 
cluster. Even closed shells are able to take advantage 
of this mechanism if they contain too few particles. If 
a small cluster approaches a strained region of a shell it 
may be able to insert itself to form a complete target 
cluster; again, any unneeded particles may be ejected. 
These events are readily visible in movies of our simula- 
tions. As an example of a typical event, a ten-particle 
closed cluster was observed to encounter a three-particle 
triangle. The triangle became loosely bound to the out- 
side of the closed shell, then over time the particles in 
the shell moved apart such that two particles from the 
triangle were able to incorporate into the shell to form 
a perfect icosahedron. The bonds with the remaining 
particle were broken, and it diffused away. 

These processes depend on a ang taking at least a mod- 
erate value for several reasons. Firstly, collisions between 
partly formed clusters are much more likely to result in 
binding interactions when the patches are wider. Sec- 
ondly, insertion of small clusters into closed shells only 
becomes possible when the patches are sufficiently wide 
that both the inserting fragment can become attached 
to the outside of the shell, despite the highly suboptimal 
angles between the respective patches, and the distortion 
required for insertion is feasible. Most importantly, the 
energy barriers for rearrangement are greatly decreased 
with wide patches, since the intermediate states can be 
stabilised by interactions between poorly aligned patches. 
These mechanisms of combination and rearrangement 
lead, for moderate values of a ang , to reasonable yields 
at long times even at very low temperatures. 

For even higher values of cr ang , yields fall off once more, 
and the kinetics in this regime are illustrated in Fig. 
|5|d). While the mechanisms of combination and rear- 
rangement remain effective (hence the number of small 
intermediates again falls off with time), the system is 
no longer so strongly driven to form correctly assembled 
clusters, since mildly misformed clusters will have similar 
energies per particle to correctly formed clusters. Excess 
particles will only rarely be ejected, and many of the 
clusters formed are oversized. 

Fig. [6] provides an alternative perspective on many of 
the effects discussed above. It shows the yields as a func- 
tion of temperature and time for each of the shapes, at 
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FIG. 5: (Colour Online) The populations of particles in clusters of different sizes as a function of number of MC cycles for 
systems of 1200 icosahedron- forming particles, averaged over 100 repetitions. The temperatures and patch widths are as 
labelled, such that (a) is at the optimal conditions for assembly, (b) is in the region of monomer starvation, (c) shows recovery 
at moderate patch widths, and (d) shows a loss of specificity at high patch widths. A log scale is used so that the very small 
populations of intermediates are discernible. States with zero population have been set to a minimum value of 0.1 to aid in 
visualisation. 



the optimal value of <r ang in each case. A number of 
effects are visible. Considering first the plot for icosahe- 
dra, Fig. |6jd) , a sharp cut-off in the early yields is seen at 
around T = 0.13, corresponding to the onset of monomer 
starvation. Yields below this cutoff recover at longer 
times, as a result of the combination and rearrange- 
ment of clusters. Similar behaviour is observed for the 
other shapes, although the sharpness of the monomer- 
starvation cutoff decreases for the smaller shapes, since 
monomer starvation affects smaller targets less seriously 
because there are a smaller number of intermediate states 
in which to get trapped. Furthermore, recovery is also 
generally easier for the smaller targets, as less rearrange- 
ment is likely to be required. The virtual absence of any 
recovery in Fig. [6je) is primarily because optimal assem- 
bly for the dodecahedron occurs at a lower value of <r ang ; 
there is some evidence of some recovery at larger values 
of cr ang in Fig. [If e). 

Fig. [7] shows the yields for each target structure as a 
function of time. The kinetics is seen to be sigmoidal, 
with an initial lag phase during which intermediates of 
increasing sizes build up in turn before the first com- 
plete clusters are formed. This is a ubiquitous feature in 



simulations of virus assembly, and indeed predicted for 
all multistep reactions.^ For larger targets the lag phase 
is longer, since the reaction has to progress through a 
larger number of intermediates. In all cases the yield 
approaches 100% at long times. 

The successful formation of dodecahedra is a partic- 
ularly notable consequence of introducing torsional con- 
straints. In the absence of torsional constraints the as- 
sembly of dodecahedra is entirely prevented by competi- 
tion with disordered aggregates, which are more thermo- 
dynamically stable at high temperatures, and which re- 
arrange only extremely slowly at lower temperatures (for 
more details see the accompanying paperi^S) However, 
since the formation of aggregates is largely prevented 
by the inclusion of torsional constraints the formation 
of dodecahedra becomes possible, and indeed approaches 
100% yield at long times. Indeed, Fig. [7] presents a pic- 
ture of the time scale for assembly increasing monoton- 
ically with target size with no obviously anomalous be- 
haviour due to particular geometric features of a target. 

We can see the effect of this lack of competition with 
aggregation in the thermodynamics of the assembly pro- 
cess. We consider first Fig. |8j which shows the heat ca- 
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FIG. 6: (Colour Online) The yields of (a) tetrahedra, (b) oc- 
tahedra, (c) cubes, (d) icosahedra and (e) dodecahedra after 
different numbers of simulation steps as a function of temper- 
ature T, in simulations of systems of 1200 particles. In each 
case the value of a ang was chosen where maximum yields were 
obtained in Fig.[l] giving a ang = 0.5 except in the case of the 
dodecahedra, where a ang = 0.35. Each data point is an av- 
erage over five simulations. All the plots use the same key, 
inset in the plot for dodecahedra. 
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FIG. 7: (Colour Online) The yields of the various target struc- 
tures as a function of time, under optimal conditions in each 
case. The temperatures used were as follows. Tetrahedra: 
T = 0.1, octahedra: T = 0.12, cubes: T = 0.09, icosahedra: 
T = 0.14, dodecahedra: T = 0.08. cr ang = 0.5 for all targets 



except the dodecahedra, where cr a 



0.35. Each line is an 



average of 100 simulations, each containing 1200 particles. 
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FIG. 8: (Colour Online) The heat capacity C v as a function 
of temperature for the formation of a single icosahedron and 
dodecahedron. The dashed lines are the analogous plots for 
systems with no torsional constraints. 



pacity C v as a function of temperature for the forma- 
tion of a single icosahedron and dodecahedron, with and 
without torsional constraints. There is a broad shoul- 
der above the main peak of the dodecahedron in the ab- 
sence of torsional constraints, because, as the tempera- 
ture is increased, the cluster first melts, before gradually 
evaporating to form a monomeric vapour. Torsional con- 
straints destabilise the disordered cluster state, leaving 
a single sharp transition between the dodecahedron and 
the monomer gas. For the icosahedron no such dramatic 
change is observed, since the icosahedron directly dis- 
assembles into monomers in both cases. Rather, there 
is just a small increase in T c corresponding to the slight 
destabilisation of the monomer gas, because it is less non- 
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FIG. 9: (Colour Online) Free energy landscapes for the forma- 
tion of a single icosahedron at a ang = 0.45 and (a) T = 0.182, 
corresponding to the peak in C v , (b) T = 0.16, the opti- 
mal temperature for assembly. The white dots indicate sub- 
clusters of the icosahedron that maximize the number of tri- 
angles for that clusters size. 



ideal. 

Two-dimensional free energy landscapes can provide a 
picture of the pathway for self-assembly. Fig. [9] shows 
that there is considerable structural order (as measured 
by the number of triangles) in the intermediates for icosa- 
hedron formation. This order is induced by torsional con- 
straints and helps the system avoid kinetic traps associ- 
ated with misformed clusters and therefore to assemble 
correctly. At the optimal conditions, the most proba- 
ble states for each cluster size are those icosahedral in- 
termediates that maximize the number of triangles (Fig. 
[9ja)). This result provides some support for those kinetic 
models that just consider the lowest-energ y p athway^ 
or a limited set of low-energy pathways! 36 * 41 1 However, 
states with fewer triangles still have significant proba- 
bilities. Furthermore, at T c the higher entropy of states 
with fewer triangles makes them more stable, and there 
is a broad valley across the free energy landscape, thus 
implying that many pathways are relevant. 

Fig. 10 'a) shows the same type of free energy landscape 
for the formation of a dodecahedron at T c . Firstly, there 



are two clear free energy minima with a barrier between 
them, confirming that the transition is now between a 
monomeric gas and the assembled dodecahedron. Sec- 
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FIG. 10: (Colour Online) Free energy landscapes for the for- 
mation of a single dodecahedron. In (a) the order parameters 
are the numbers of particles and pentagons in the largest clus- 
ter, while in (b) they are the potential energy and a measure 
of radial disorder, described in the text. The insets show 
the analogous plots for systems with no torsional constraints, 
where the axes of the insets span the same parameter ranges 
as the main plots. The white dots in (a) indicate sub-clusters 
of the dodecahedron that maximize the number of pentagons 
for that clusters size. 



ondly, similar to the icosahedron, there is considerable 
structural order along the pathway for assembly, and so 
as clusters grow they are 'guided' towards the dodecahe- 
dral target. This free energy landscape contrasts sharply 
with the analogous landscape for a system without tor- 
sional constraints (see inset), for which disordered aggre- 
gates are also stable and block the formation of the target 
structure. 

A different view of the free energy landscape is pro- 
vided in Fig. 10 'b) where the order parameters are the 
potential energy and a "radial order parameter" that is 
defined as the standard deviation in the distance of par- 
ticles from the centre of mass, and takes a value of zero 
when all particles lie on the surface of a sphere (as in a 
perfect dodecahedron). In the absence of torsional con- 
straints the system is able to access very disordered states 
with little energy penalty, and indeed can reach low en- 
ergies while remaining very far from the dodecahedral 
structure. Once torsional constraints are included the 
landscape takes on a classic funnel shapep^ such that as 
the energy of the system decreases it is directed into the 
dodecahedral structure. 
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IV. CONCLUSIONS 

We have extended a patchy particl e mod el that can 
self-assemble into monodisperse clusters^^to include a 
torsional component in the interparticle potential, better 
mimicking the nature of protein-protein interactions. We 
found that similar principles operate in both these mod- 
els. Firstly, in both cases there is a temperature window 
over which successful assembly can occur because at too 
high temperatures the target structure is thermodynam- 
ically unstable and at low temperatures the system be- 
comes kinetically trapped. Similarly, there is a limited 
range of patch specificities associated with successful as- 
sembly, because when the patches are too wide alterna- 
tive configurations become possible, and when the patch- 
patch interactions are too specific growth is too slow. 

However, there are also clear differences between the 
models. In particular, the torsional constraints greatly 
reduce the set of erroneous structures that can be formed, 
and thus remove the competition with large disordered 
aggregates that was a feature of the original model. So 
now, rather than large disordered aggregates, the alter- 
native configurations at larger a are malformed shell-like 
clusters, and at lower temperatures the system becomes 
trapped in incomplete clusters because of monomer star- 
vation due to the rapid rate of cluster nucleation. Both 
these features make the behaviour of the model much 
more similar to that for virus capsids, as would be ex- 
pected from the more protein-like nature of the interac- 
tions. 

This lack of competition from aggregates enables the 
assembly of larger target clusters than was possible in the 
original model. In particular, the dodecahedron which 
proved impossible to assemble without torsional inter- 
actions readily assembles for the current model. The 
free energy landscapes for the respective systems dra- 
matically illustrate these differences. With torsional con- 
straints, the free energy landscapes now have a funnel-like 
topography guiding the system towards the target clus- 
ter, and intermediate-sized clusters exhibit a significant 



amount of the structural order that is characteristic of 
the target, i.e. as the clusters grow, they grow with the 
correct structure. 

Interestingly, the assembly behaviour for the different 
targets is very similar, save for the general effects of tar- 
get size, e.g. the time taken to achieve assembly increases 
monotonically with size, and the effects of monomer 
starvation become more pronounced for larger targets. 
This similarity contrasts with the original non-torsional 
model, where the dependence of the propensity to ag- 
gregate on the patch geometry leads to large differences 
in behaviour for the different targets. Notably, the abil- 
ity of all the targets to readily assemble for the current 
model suggests that proteins are less limited in the geo- 
metric forms into which they can assemble, as also seems 
apparent from the diverse range of biological structures 
and machines made up of proteins. 

If synthetic particles are to begin to mimic this bi- 
ological repertoire, our results suggest that torsionally- 
specific interactions between the particles would be re- 
quired. However, it is not apparent how such orienta- 
tional specificity might be achieved for patchy nanopar- 
ticles and colloids. By contrast, the DNA units con- 
structed by Mao and coworkers have both 'valency' and 
control over the relative orientation of the units. Thus, 
our results help to understand why tetrahedra, cubes, do- 
decahedra and truncated icosahedra can be successfully 
assembled from 3-armed units pES anc [ icosahedra from 
5-armed units. 59 Although far from spherical, the cur- 
rent model may not provide such a bad representation 
for these DNA systems, especially if internal degrees of 
freedom for the particles could be added to mimic the 
flexibility of the arms. 
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